##########################################################
### Disorganized Political Violence:                   ###
### A Demonstration Case of Temperature and Insurgency ###
### Replication Code:                                  ###
### Temperature Deviation Results:                     ###
### Mean Temperatures Robustness Check                 ###
### Andrew Shaver, Alex Bollfrass                      ###
### Date: 07/04/2022                                   ###
##########################################################

# Note: these results were generated using:
# 1) MacBook Pro (16-inch, 2021); Chip: Apple M1 Max; Memory: 64 GB; running
# macOS Monterey.
# 2) R 4.1.2 GUI 1.77 High Sierra build (8007)
# 3) RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS
# Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
# (Please see below for specific package versions.)

# Load required packages:
library(MASS)
library(lubridate)
library(DataCombine)
library(broom)
library(data.table)
library(radiant.model)
library(lfe)
library(ggplot2)
library(ggpubr)
library(stargazer)
library(arm)

packageVersion("MASS")
# ‘7.3.55’
packageVersion("lubridate")
# ‘1.8.0’
packageVersion("DataCombine")
# ‘0.2.21’
packageVersion("broom")
# ‘0.7.12’
packageVersion("data.table")
# ‘1.14.2’
packageVersion("radiant.model")
# ‘1.4.4’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("ggplot2")
# ‘3.3.6’
packageVersion("ggpubr")
# ‘0.4.0’
packageVersion("stargazer")
# ‘5.2.3’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("arm")
# ‘1.12.2’

# Load data:
# For primary temperature-violence analysis:
# Daily analysis:
iq_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/iraq_baghdad_basrah_processed_data.rds")
af_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/afghanstan_viol_districts_processed_data.rds")
# For attitude analysis:
bg_sur <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/baghdad_survey.rds")
# For median analysis:
iq_median_temp <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/Affective_Motivations/MajorRevision/R3_Submission/ProcessedData/iq_daily_median.rds")

# Next, construct temperature deviations measure:
# Write function:
deviation <- function(x){
  admin_list <- list()
  library(locfit)
  library(data.table)
  for(i in 1:length(unique(panel[, x]))){
    temp <- panel[panel[, x] %in% unique(panel[, x])[i], ]
    temp$nday <- 1:nrow(temp)
    t <- locfit(temp ~ lp(nday, nn = .05), data=temp)
    temp$temp_hat <- predict(t, newdata = temp$nday)
    temp$temp_dev <- temp$temp - temp$temp_hat
    admin_list[[i]] <- temp
    rm(temp)
  }
  admin_list <- rbindlist(admin_list)
  return(admin_list)
}
panel <- iq_sig1
iq_sig1 <- deviation("governorate")

panel <- af_sig1
af_sig1 <- deviation("station")

iq_sig1 <- as.data.frame(iq_sig1)
af_sig1 <- as.data.frame(af_sig1)

# For temperature deviations for the Baghdad survey data, given that the unit is not the
# match deviations from panel:
# Columns to keep:
to_keep <- c("date", "temp_dev", "temp_1", "temp_2", "temp_3", "temp_4", 
             "temp_5", "temp_6", "temp_7", "daylight")
bg_sur <- plyr::join(bg_sur, iq_sig1[iq_sig1$governorate %in% "Baghdad", to_keep], by = "date")

# Create survey dataset removing top income earners:
bg_sur.s3.l <- bg_sur[bg_sur$income.w < 75000, ]

# Next, generate temperature deviations for each context:
# Set temperature range for all:
temperature_range <- 65:95

###
### IRAQ.
###

# Baghdad and Basrah:
# Unconstrained Violence:
# OLS model results:

# Check the normality of the errors:
iq_sig1$threshold <- ifelse(iq_sig1$temp > 75, 1, 0)
test1 <- lm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(governorate), data = iq_sig1)
test2 <- lm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(governorate), data = iq_sig1)

test1.stdres = rstandard(test1)
qqnorm(test1.stdres)
qqline(test1.stdres)

test2.stdres = rstandard(test2)
qqnorm(test2.stdres)
qqline(test2.stdres)

temperature_deviation_iq_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                                  temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained |
                                  week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(unconstrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols <- temperature_deviation_iq_ols(iq_sig1, 0.95)
iq_90_ols <- temperature_deviation_iq_ols(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols$change <- (sd(iq_sig1$temp_dev, na.rm = T) * iq_95_ols$estimate) / mean(log(iq_sig1$unconstrained+1))

# For plotting Baghdad-Basrah results:
temp_ranges <- unique(iq_95_ols$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges)*2-2), by=2)){
  temp_ranges <- append(temp_ranges, " ", after=i)
}

# Plot:
p_iq_ols <- ggplot() + 
  
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_ols[iq_95_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_ols[iq_90_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_ols[iq_90_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 65:95, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# Negative binomial model results:

# Calculate magnitude in terms of percentage deviation from the mean:
new_data_iq <- data.frame(temp_dev = c(0, sd(iq_sig1$temp_dev)),
                          threshold = 0,
                          
                          prcp = mean(iq_sig1$prcp, na.rm = T), 
                          wdsp = mean(iq_sig1$wdsp, na.rm = T), 
                          dewp = mean(iq_sig1$dewp, na.rm = T), 
                          visib = mean(iq_sig1$visib, na.rm = T), 
                          mxspd = mean(iq_sig1$mxspd, na.rm = T), 
                          daylight = mean(iq_sig1$daylight, na.rm = T), 
                          hoursofpower = mean(iq_sig1$hoursofpower, na.rm = T), 
                          
                          idf = mean(iq_sig1$idf), 
                          ied = mean(iq_sig1$ied),  
                          constrained = mean(iq_sig1$constrained),  
                          unconstrained = mean(iq_sig1$unconstrained),  
                          
                          temp_1 = mean(iq_sig1$temp_1, na.rm = T),
                          temp_2 = mean(iq_sig1$temp_2, na.rm = T),
                          temp_3 = mean(iq_sig1$temp_3, na.rm = T),
                          temp_4 = mean(iq_sig1$temp_4, na.rm = T),
                          temp_5 = mean(iq_sig1$temp_5, na.rm = T),
                          temp_6 = mean(iq_sig1$temp_6, na.rm = T),
                          temp_7 = mean(iq_sig1$temp_7, na.rm = T),
                          
                          unconstrained_1 = mean(iq_sig1$unconstrained_1, na.rm = T),
                          unconstrained_2 = mean(iq_sig1$unconstrained_2, na.rm = T),
                          unconstrained_3 = mean(iq_sig1$unconstrained_3, na.rm = T),
                          unconstrained_4 = mean(iq_sig1$unconstrained_4, na.rm = T),
                          unconstrained_5 = mean(iq_sig1$unconstrained_5, na.rm = T),
                          unconstrained_6 = mean(iq_sig1$unconstrained_6, na.rm = T),
                          unconstrained_7 = mean(iq_sig1$unconstrained_7, na.rm = T),
                          
                          constrained_1 = mean(iq_sig1$constrained_1, na.rm = T),
                          constrained_2 = mean(iq_sig1$constrained_2, na.rm = T),
                          constrained_3 = mean(iq_sig1$constrained_3, na.rm = T),
                          constrained_4 = mean(iq_sig1$constrained_4, na.rm = T),
                          constrained_5 = mean(iq_sig1$constrained_5, na.rm = T),
                          constrained_6 = mean(iq_sig1$constrained_6, na.rm = T),
                          constrained_7 = mean(iq_sig1$constrained_7, na.rm = T),
                          
                          idf_1 = mean(iq_sig1$idf_1, na.rm = T),
                          idf_2 = mean(iq_sig1$idf_2, na.rm = T),
                          idf_3 = mean(iq_sig1$idf_3, na.rm = T),
                          idf_4 = mean(iq_sig1$idf_4, na.rm = T),
                          idf_5 = mean(iq_sig1$idf_5, na.rm = T),
                          idf_6 = mean(iq_sig1$idf_6, na.rm = T),
                          idf_7 = mean(iq_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(iq_sig1$ied_1, na.rm = T),
                          ied_2 = mean(iq_sig1$ied_2, na.rm = T),
                          ied_3 = mean(iq_sig1$ied_3, na.rm = T),
                          ied_4 = mean(iq_sig1$ied_4, na.rm = T),
                          ied_5 = mean(iq_sig1$ied_5, na.rm = T),
                          ied_6 = mean(iq_sig1$ied_6, na.rm = T),
                          ied_7 = mean(iq_sig1$ied_7, na.rm = T),
                          
                          week = median(iq_sig1$week),
                          month = median(iq_sig1$month),
                          governorate = "Baghdad")

temperature_deviation_iq_nb <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                                 temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                                 as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$unconstrained)
    
    reg_temp_month.nb <- glm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
             as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(unconstrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                                   as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"

  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"

  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"

  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb <- temperature_deviation_iq_nb(iq_sig1, 0.95)
iq_90_nb <- temperature_deviation_iq_nb(iq_sig1, 0.90)

# Plot:

p_iq_nb <- ggplot() + 
  
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = iq_95_nb[iq_95_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = iq_90_nb[iq_90_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = iq_90_nb[iq_90_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "blue4") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = iq_95_nb[iq_95_nb$time %in% "week", ]$gamma, labels=scales::percent(iq_95_nb[iq_95_nb$time %in% "week", ]$change, accuracy = 0.01))) +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 65:95, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# Constrained Violence:
# OLS model results:
temperature_deviation_iq_ols_con <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(constrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(constrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(constrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(constrained+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols_con <- temperature_deviation_iq_ols_con(iq_sig1, 0.95)
iq_90_ols_con <- temperature_deviation_iq_ols_con(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols_con$change <- (sd(iq_sig1$temp_dev, na.rm = T) * iq_95_ols_con$estimate) / mean(log(iq_sig1$constrained+1))

# Plot:
p_iq_ols_con <- ggplot() + 
  
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "green4") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") + 
  geom_errorbar(data = iq_95_ols_con[iq_95_ols_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "green4", linetype = "dashed") +
  geom_point(data = iq_90_ols_con[iq_90_ols_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") +
  geom_errorbar(data = iq_90_ols_con[iq_90_ols_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "green4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 65:95, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

temperature_deviation_iq_nb_con <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(constrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$constrained)
    
    reg_temp_month.nb <- glm(constrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(constrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                                   as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(constrained ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                                    as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb_con <- temperature_deviation_iq_nb_con(iq_sig1, 0.95)
iq_90_nb_con <- temperature_deviation_iq_nb_con(iq_sig1, 0.90)

# Plot:

p_iq_nb_con <- ggplot() + 
  
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "green4") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") + 
  geom_errorbar(data = iq_95_nb_con[iq_95_nb_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "green4", linetype = "dashed") +
  geom_point(data = iq_90_nb_con[iq_90_nb_con$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "green4") +
  geom_errorbar(data = iq_90_nb_con[iq_90_nb_con$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "green4") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = iq_95_nb_con[iq_95_nb_con$time %in% "week", ]$gamma, labels=scales::percent(iq_95_nb_con[iq_95_nb_con$time %in% "week", ]$change, accuracy = 0.01))) + 
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 65:95, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

# ied Attacks:
# OLS model results:
temperature_deviation_iq_ols_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(log(ied+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + governorate | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(ied+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                             unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + governorate | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(ied+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained |
                                 week + governorate | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(ied+1) ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained |
                                  month + governorate | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_ols_ied <- temperature_deviation_iq_ols_ied(iq_sig1, 0.95)
iq_90_ols_ied <- temperature_deviation_iq_ols_ied(iq_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
iq_95_ols_ied$change <- (sd(iq_sig1$temp_dev, na.rm = T) * iq_95_ols_ied$estimate) / mean(log(iq_sig1$ied+1))

# Plot:
p_iq_ols_ied <- ggplot() + 
  
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = iq_95_ols_ied[iq_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = iq_90_ols_ied[iq_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "firebrick") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.02, 0.0275) + 
  annotate(geom = "text", x = 65:95, y = -0.0165,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

temperature_deviation_iq_nb_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.nb <- glm(ied ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                              unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint_robust(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint_robust(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_iq[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_iq[1,], type = "response")) /  mean(iq_sig1$ied)
    
    reg_temp_month.nb <- glm(ied ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                               unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint_robust(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint_robust(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(ied ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                                   as.factor(week) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint_robust(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint_robust(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(ied ~ temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                                    as.factor(month) + as.factor(governorate), data = temp, family = quasipoisson)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint_robust(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint_robust(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

iq_95_nb_ied <- temperature_deviation_iq_nb_ied(iq_sig1, 0.95)
iq_90_nb_ied <- temperature_deviation_iq_nb_ied(iq_sig1, 0.90)

# Plot:

p_iq_nb_ied <- ggplot() + 
  
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = iq_90_nb_ied[iq_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "firebrick") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ]$gamma, labels=scales::percent(iq_95_nb_ied[iq_95_nb_ied$time %in% "week", ]$change, accuracy = 0.01))) + 
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.035, 0.05) + 
  annotate(geom = "text", x = 65:95, y = -0.03,
           label = temp_ranges,
           size = 3.5, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=12),
        axis.text.x.top = element_text(angle = 90, hjust = 0, size=10),
        axis.title=element_text(size=12,face="bold"))

iraq_combined <- ggarrange(p_iq_nb + xlab("") + ylab("Quasi-Poisson"), 
          p_iq_nb_con + xlab("") + ylab(""), 
          p_iq_nb_ied + xlab("") + ylab(""),  
          p_iq_ols + xlab("") + ylab("OLS"),
          p_iq_ols_con + xlab("") + ylab(""), 
          p_iq_ols_ied + xlab("") + ylab(""), 
          ncol = 3, nrow=2, labels  = "AUTO")

# PRODUCES FIGURE 12 IN PAPER:
annotate_figure(iraq_combined,
                top = text_grob(expression(rho), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*" (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### AFGHANISTAN.
###

# Afghanistan districts:
# OLS:

# Check the normality of the errors:
af_sig1$threshold <- ifelse(af_sig1$temp > 75, 1, 0)
test1 <- lm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                        temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                        df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                        idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                        ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                        as.factor(week) + as.factor(station), data = af_sig1)
test2 <- lm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
              as.factor(week) + as.factor(station), data = af_sig1)

test1.stdres = rstandard(test1)
qqnorm(test1.stdres)
qqline(test1.stdres)

test2.stdres = rstandard(test2)
qqnorm(test2.stdres)
qqline(test2.stdres)

# DF:
# Log:
temperature_deviation_af_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > (i + (min(temperature_range)-1)), 1, 0)
    
    reg_temp_week <- felm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                                 temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                                 df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                                 idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                                 ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                                 week + station | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
  
    reg_temp_month <- felm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + station | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))

    reg_temp_week_pars <- felm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied |
                                 week + station | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))

    reg_temp_month_pars <- felm(log(df+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied |
                                  month + station | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"

  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"

  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"

  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_ols <- temperature_deviation_af_ols(af_sig1, 0.95)
af_90_ols <- temperature_deviation_af_ols(af_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
af_95_ols$change <- (sd(af_sig1$temp_dev, na.rm = T) * af_95_ols$estimate) / mean(log(af_sig1$df+1))

# For plotting Afghan district results:
temp_ranges_af <- unique(af_95_ols$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges_af)*2-2), by=2)){
  temp_ranges_af <- append(temp_ranges_af, " ", after=i)
}

p_af_ols <- ggplot() + 
  
  geom_line(data = af_95_ols[af_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = af_95_ols[af_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_ols[af_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_ols[af_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_ols[af_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols[af_95_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols[af_90_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols[af_90_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols[af_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols[af_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols[af_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols[af_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols[af_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols[af_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols[af_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols[af_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols[af_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = af_95_ols[af_95_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = af_90_ols[af_90_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = af_90_ols[af_90_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.01, 0.0165) + 
  annotate(geom = "text", x = 65:95, y = -0.0085,
           label = temp_ranges_af,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# Negative binomial model:

# Calculate magnitude in terms of percentage deviation from the mean:
new_data_af <- data.frame(temp_dev = c(0, sd(af_sig1$temp_dev)),
                          threshold = 0,
                        
                          prcp = mean(af_sig1$prcp, na.rm = T), 
                          wdsp = mean(af_sig1$wdsp, na.rm = T), 
                          dewp = mean(af_sig1$dewp, na.rm = T), 
                          visib = mean(af_sig1$visib, na.rm = T), 
                          mxspd = mean(af_sig1$mxspd, na.rm = T), 
                          daylight = mean(af_sig1$daylight, na.rm = T), 

                          idf = mean(af_sig1$idf), 
                          ied = mean(af_sig1$ied),  
                          df = mean(af_sig1$df), 
                          
                          temp_1 = mean(af_sig1$temp_1, na.rm = T),
                          temp_2 = mean(af_sig1$temp_2, na.rm = T),
                          temp_3 = mean(af_sig1$temp_3, na.rm = T),
                          temp_4 = mean(af_sig1$temp_4, na.rm = T),
                          temp_5 = mean(af_sig1$temp_5, na.rm = T),
                          temp_6 = mean(af_sig1$temp_6, na.rm = T),
                          temp_7 = mean(af_sig1$temp_7, na.rm = T),
                          
                          df_1 = mean(af_sig1$df_1, na.rm = T),
                          df_2 = mean(af_sig1$df_2, na.rm = T),
                          df_3 = mean(af_sig1$df_3, na.rm = T),
                          df_4 = mean(af_sig1$df_4, na.rm = T),
                          df_5 = mean(af_sig1$df_5, na.rm = T),
                          df_6 = mean(af_sig1$df_6, na.rm = T),
                          df_7 = mean(af_sig1$df_7, na.rm = T),
                          
                          idf_1 = mean(af_sig1$idf_1, na.rm = T),
                          idf_2 = mean(af_sig1$idf_2, na.rm = T),
                          idf_3 = mean(af_sig1$idf_3, na.rm = T),
                          idf_4 = mean(af_sig1$idf_4, na.rm = T),
                          idf_5 = mean(af_sig1$idf_5, na.rm = T),
                          idf_6 = mean(af_sig1$idf_6, na.rm = T),
                          idf_7 = mean(af_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(af_sig1$ied_1, na.rm = T),
                          ied_2 = mean(af_sig1$ied_2, na.rm = T),
                          ied_3 = mean(af_sig1$ied_3, na.rm = T),
                          ied_4 = mean(af_sig1$ied_4, na.rm = T),
                          ied_5 = mean(af_sig1$ied_5, na.rm = T),
                          ied_6 = mean(af_sig1$ied_6, na.rm = T),
                          ied_7 = mean(af_sig1$ied_7, na.rm = T),
                          
                          week = median(af_sig1$week),
                          month = median(af_sig1$month),
                          station = "409900")

temperature_deviation_af_nb <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > (i + (min(temperature_range)-1)), 1, 0)
    
    reg_temp_week.nb <- glm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                            as.factor(week) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint.default(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint.default(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_af[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_af[1,], type = "response")) /  mean(af_sig1$df)
    
    reg_temp_month.nb <- glm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                             as.factor(month) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint.default(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint.default(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                                 as.factor(week) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint.default(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint.default(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(df ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                                  as.factor(month) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint.default(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint.default(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_nb <- temperature_deviation_af_nb(af_sig1, 0.95)
af_90_nb <- temperature_deviation_af_nb(af_sig1, 0.90)

p_af_qp <- ggplot() + 
  
  geom_line(data = af_95_nb[af_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = af_95_nb[af_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_nb[af_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_nb[af_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_nb[af_95_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb[af_95_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb[af_90_nb$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb[af_90_nb$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb[af_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb[af_95_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb[af_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb[af_90_nb$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb[af_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb[af_95_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb[af_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb[af_90_nb$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb[af_95_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") + 
  geom_errorbar(data = af_95_nb[af_95_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = af_90_nb[af_90_nb$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "blue4") +
  geom_errorbar(data = af_90_nb[af_90_nb$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "blue4") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = af_95_nb[af_95_nb$time %in% "week", ]$gamma, labels=scales::percent(af_95_nb[af_95_nb$time %in% "week", ]$change, accuracy = 0.01))) + 
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.01, 0.0165) + 
  annotate(geom = "text", x = 65:95, y = -0.0085,
           label = temp_ranges_af,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# IED:
temperature_deviation_af_ols_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > (i + (min(temperature_range)-1)), 1, 0)
    
    reg_temp_week <- felm(log(ied+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                            idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                            week + station | 0 | 0, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(log(ied+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                             df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                             idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                             ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                             month + station | 0 | 0, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(log(ied+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df |
                                 week + station | 0 | 0, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(log(ied+1) ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df |
                                  month + station | 0 | 0, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, se.type = "robust", conf.int = T, conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_ols_ied <- temperature_deviation_af_ols_ied(af_sig1, 0.95)
af_90_ols_ied <- temperature_deviation_af_ols_ied(af_sig1, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
af_95_ols_ied$change <- (sd(af_sig1$temp_dev, na.rm = T) * af_95_ols_ied$estimate) / mean(log(af_sig1$ied+1))

p_af_ols_ied <- ggplot() + 
  
  geom_line(data = af_95_ols_ied[af_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = af_95_ols_ied[af_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_ols_ied[af_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_ols_ied[af_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_ols_ied[af_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols_ied[af_95_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols_ied[af_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols_ied[af_90_ols_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols_ied[af_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols_ied[af_95_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols_ied[af_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols_ied[af_90_ols_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols_ied[af_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_ols_ied[af_95_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_ols_ied[af_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_ols_ied[af_90_ols_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_ols_ied[af_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = af_95_ols_ied[af_95_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = af_90_ols_ied[af_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = af_90_ols_ied[af_90_ols_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "firebrick") +
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.01, 0.0165) + 
  annotate(geom = "text", x = 65:95, y = -0.0085,
           label = temp_ranges_af,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# QP:
temperature_deviation_af_nb_ied <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > (i + (min(temperature_range)-1)), 1, 0)
    
    reg_temp_week.nb <- glm(ied ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                              temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                              df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                              idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                              ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                              as.factor(week) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint.default(reg_temp_week.nb, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint.default(reg_temp_week.nb, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.nb, new_data_af[2,], type = "response") - 
                                    predict(reg_temp_week.nb, new_data_af[1,], type = "response")) /  mean(af_sig1$df)
    
    reg_temp_month.nb <- glm(ied ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                               temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                               idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                               ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                               as.factor(month) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint.default(reg_temp_month.nb, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint.default(reg_temp_month.nb, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_week_pars.nb <- glm(ied ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                                   as.factor(week) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint.default(reg_temp_week_pars.nb, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint.default(reg_temp_week_pars.nb, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month_pars.nb <- glm(ied ~ temp_dev + temp_dev*threshold + prcp + wdsp + dewp + visib + mxspd + daylight + idf + df +
                                    as.factor(month) + as.factor(station), family = quasipoisson, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.nb) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint.default(reg_temp_month_pars.nb, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint.default(reg_temp_month_pars.nb, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

af_95_nb_ied <- temperature_deviation_af_nb_ied(af_sig1, 0.95)
af_90_nb_ied <- temperature_deviation_af_nb_ied(af_sig1, 0.90)

p_af_qp_ied <- ggplot() + 
  
  geom_line(data = af_95_nb_ied[af_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "firebrick") +
  geom_line(data = af_95_nb_ied[af_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_nb_ied[af_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = af_95_nb_ied[af_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = af_95_nb_ied[af_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb_ied[af_95_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb_ied[af_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb_ied[af_90_nb_ied$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb_ied[af_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb_ied[af_95_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb_ied[af_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb_ied[af_90_nb_ied$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb_ied[af_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_95_nb_ied[af_95_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = af_90_nb_ied[af_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = af_90_nb_ied[af_90_nb_ied$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = af_95_nb_ied[af_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") + 
  geom_errorbar(data = af_95_nb_ied[af_95_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "firebrick", linetype = "dashed") +
  geom_point(data = af_90_nb_ied[af_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, color = "firebrick") +
  geom_errorbar(data = af_90_nb_ied[af_90_nb_ied$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "firebrick") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = af_95_nb_ied[af_95_nb_ied$time %in% "week", ]$gamma, labels=scales::percent(af_95_nb_ied[af_95_nb_ied$time %in% "week", ]$change, accuracy = 0.01))) + 
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.01, 0.0165) + 
  annotate(geom = "text", x = 65:95, y = -0.0085,
           label = temp_ranges_af,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

afghanistan_combined <- ggarrange(p_af_qp + xlab("") + ylab("Quasi-Poisson"),
          p_af_qp_ied + xlab("") + ylab(""),
          p_af_ols + xlab("") + ylab("OLS"),
          p_af_ols_ied + xlab("") + ylab(""),
          ncol = 2, nrow=2, labels  = "AUTO")

# PRODUCES FIGURE 13 IN PAPER:
annotate_figure(afghanistan_combined,
                top = text_grob(expression(rho), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*" (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### BAGHDAD SURVEYS.
###

# OLS model results:
temperature_deviation_bg_ols <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week <- felm(vs_mnf ~ temp_dev*threshold + age + work_wkhr +income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                            p3_electr_binary +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                            education + mahala + date_month | 0 | mahala, data = temp)

    reg_list_week[[i]] <- tidy(reg_temp_week, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month <- felm(vs_mnf ~ temp_dev*threshold + age + work_wkhr +income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                             p3_electr_binary + 
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                             education + block + date_month | 0 | block, data = temp)
    
    reg_list_month[[i]] <- tidy(reg_temp_month, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_week_pars <- felm(vs_mnf ~ temp_dev*threshold + age + work_wkhr +income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight |
                                 p3_electr_binary +
                                 education + mahala + date_month | 0 | mahala, data = temp)

    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
    reg_temp_month_pars <- felm(vs_mnf ~ temp_dev*threshold + age + work_wkhr +income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                                p3_electr_binary |
                                education + block + date_month | 0 | block, data = temp)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars, conf.int = T, se.type = "robust", conf.level = z) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    
  }
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"

  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

bg_95_ols <- temperature_deviation_bg_ols(bg_sur.s3.l, 0.95)
bg_90_ols <- temperature_deviation_bg_ols(bg_sur.s3.l, 0.90)

# Calculate magnitude in terms of percentage deviation from the mean:
bg_95_ols$change <- (sd(bg_sur.s3.l$temp_dev, na.rm = T) * bg_95_ols$estimate) / mean(bg_sur.s3.l$vs_mnf)

# For plotting Baghdad survey results:
temp_ranges_bg <- unique(bg_95_ols$bin)[c(TRUE, FALSE)]
for(i in seq(1,(length(temp_ranges_bg)*2-2), by=2)){
  temp_ranges_bg <- append(temp_ranges_bg, " ", after=i)
}

p_bg_ols <- ggplot() +
  
  geom_line(data = bg_95_ols[bg_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = bg_95_ols[bg_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_line(data = bg_95_ols[bg_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size =  0.5, alpha = .35, color = "gray60") +
  geom_line(data = bg_95_ols[bg_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +

  geom_point(data = bg_95_ols[bg_95_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, alpha = 1, color = "blue4") +
  geom_errorbar(data = bg_95_ols[bg_95_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = bg_90_ols[bg_90_ols$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, alpha = 1, color = "blue4") +
  geom_errorbar(data = bg_90_ols[bg_90_ols$time %in% "week", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = 1, size = 1, color = "blue4") +
  
  geom_point(data = bg_95_ols[bg_95_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_ols[bg_95_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_ols[bg_90_ols$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_ols[bg_90_ols$time %in% "month", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 1, color = "gray60") +

  geom_point(data = bg_95_ols[bg_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_ols[bg_95_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_ols[bg_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_ols[bg_90_ols$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = bg_95_ols[bg_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_ols[bg_95_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_ols[bg_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_ols[bg_90_ols$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = conf.low, ymax = conf.high), alpha = .35, size = .75, color = "gray60") +

  ylab(expression(beta~Delta*"Temperature")) +
  xlab(expression(gamma*" (°F)")) +
  ylim(-0.004, 0.015) +
  annotate(geom = "text", x = 65:95, y = -0.00325,
           label = temp_ranges_bg,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") +
  theme_bw() +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# Logit model results:

# Calculate magnitude in terms of percentage deviation from the mean:
new_data_bg <- data.frame(temp_dev = c(0, sd(bg_sur.s3.l$temp_dev)),
                          threshold = 0,
                          
                          prcp = mean(bg_sur.s3.l$prcp, na.rm = T), 
                          wdsp = mean(bg_sur.s3.l$wdsp, na.rm = T), 
                          dewp = mean(bg_sur.s3.l$dewp, na.rm = T), 
                          visib = mean(bg_sur.s3.l$visib, na.rm = T), 
                          mxspd = mean(bg_sur.s3.l$mxspd, na.rm = T), 
                          daylight = mean(bg_sur.s3.l$daylight, na.rm = T), 
                          
                          age = mean(bg_sur.s3.l$age, na.rm = T),  
                          work_wkhr = mean(bg_sur.s3.l$work_wkhr),  
                          income.w = mean(bg_sur.s3.l$income.w),  
                          hhld_size = mean(bg_sur.s3.l$hhld_size), 
                          m_sig1_7 = mean(bg_sur.s3.l$m_sig1_7),  
                          p3_electr_binary = mean(bg_sur.s3.l$p3_electr_binary, na.rm = T),
                          
                          temp_1 = mean(bg_sur.s3.l$temp_1, na.rm = T),
                          temp_2 = mean(bg_sur.s3.l$temp_2, na.rm = T),
                          temp_3 = mean(bg_sur.s3.l$temp_3, na.rm = T),
                          temp_4 = mean(bg_sur.s3.l$temp_4, na.rm = T),
                          temp_5 = mean(bg_sur.s3.l$temp_5, na.rm = T),
                          temp_6 = mean(bg_sur.s3.l$temp_6, na.rm = T),
                          temp_7 = mean(bg_sur.s3.l$temp_7, na.rm = T),
                          
                          education = "Baccalaureate",
                          date_month = median(bg_sur.s3.l$date_month),
                          date_week = median(bg_sur.s3.l$date_week),
                          mahala = "Kadhimyah",
                          block = "1001")

temperature_deviation_bg_logit <- function(x, z){ 
  # x == dataset
  # z == confidence level
  reg_list_week <- list()
  reg_list_month <- list()
  reg_list_week_pars <- list()
  reg_list_month_pars <- list()
  
  for(i in 1:length(temperature_range)){
    temp <- x
    temp$threshold <- ifelse(temp$temp > i + (min(temperature_range)-1), 1, 0)
    
    reg_temp_week.l <- bayesglm(vs_mnf ~ temp_dev*threshold + age + work_wkhr +income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + 
                                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                                   p3_electr_binary + 
                                   as.factor(education) + as.factor(mahala) + as.factor(date_month), 
                                 data = temp, family = binomial(link = "logit"), drop.unused.levels = F)
    
    reg_list_week[[i]] <- tidy(reg_temp_week.l) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week[[i]]$lower <- confint.default(reg_temp_week.l, level = z)[2,1]
    reg_list_week[[i]]$upper <- confint.default(reg_temp_week.l, level = z)[2,2]
    reg_list_week[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week[[i]]$change <- (predict(reg_temp_week.l, new_data_bg[2,], type = "response") - 
                                     predict(reg_temp_week.l, new_data_bg[1,], type = "response")) /  mean(bg_sur.s3.l$vs_mnf)
    
    reg_temp_week_pars.l <- bayesglm(vs_mnf ~ temp_dev*threshold + age + work_wkhr +income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd +
                                        p3_electr_binary + 
                                        as.factor(education) + as.factor(mahala) + as.factor(date_month), 
                                      data = temp, family = binomial(link = "logit"), drop.unused.levels = F)
    
    reg_list_week_pars[[i]] <- tidy(reg_temp_week_pars.l) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_week_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_week_pars[[i]]$lower <- confint.default(reg_temp_week_pars.l, level = z)[2,1]
    reg_list_week_pars[[i]]$upper <- confint.default(reg_temp_week_pars.l, level = z)[2,2]
    reg_list_week_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_week_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_week_pars[[i]]$change <- NA
    
    reg_temp_month.l <- bayesglm(vs_mnf ~ temp_dev*threshold + age + work_wkhr +income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + 
                             temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                               p3_electr_binary + 
                               as.factor(education) + as.factor(block) + as.factor(date_month), 
                             data = temp, family = binomial(link = "logit"), drop.unused.levels = F)
    
    reg_list_month[[i]] <- tidy(reg_temp_month.l) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month[[i]]$lower <- confint.default(reg_temp_month.l, level = z)[2,1]
    reg_list_month[[i]]$upper <- confint.default(reg_temp_month.l, level = z)[2,2]
    reg_list_month[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month[[i]]$change <- NA
    
    reg_temp_month_pars.l <- bayesglm(vs_mnf ~ temp_dev*threshold + age + work_wkhr +income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd +
                                        p3_electr_binary + 
                                        as.factor(education) + as.factor(block) + as.factor(date_month), 
                                      data = temp, family = binomial(link = "logit"), drop.unused.levels = F)
    
    reg_list_month_pars[[i]] <- tidy(reg_temp_month_pars.l) %>% filter(term %in% "temp_dev") %>% mutate(model = paste("Model",i,sep=" "))
    reg_list_month_pars[[i]]$bin <- paste("[", round(min(x$temp)), ", ", as.character(i + (min(temperature_range)-1)), ")", sep="")
    reg_list_month_pars[[i]]$lower <- confint.default(reg_temp_month_pars.l, level = z)[2,1]
    reg_list_month_pars[[i]]$upper <- confint.default(reg_temp_month_pars.l, level = z)[2,2]
    reg_list_month_pars[[i]]$gamma <- i + (min(temperature_range)-1)
    reg_list_month_pars[[i]]$temperature_pct <- ecdf(x$temp)(i + (min(temperature_range)-1))
    reg_list_month_pars[[i]]$change <- NA
    
  }
  
  week_fe <- rbindlist(reg_list_week)
  week_fe$time <- "week"
  
  week_fe_pars <- rbindlist(reg_list_week_pars)
  week_fe_pars$time <- "week_pars"
  
  month_fe <- rbindlist(reg_list_month)
  month_fe$time <- "month"
  
  month_fe_pars <- rbindlist(reg_list_month_pars)
  month_fe_pars$time <- "month_pars"
  
  week_month_fe <- rbind(week_fe, week_fe_pars, month_fe, month_fe_pars)
  return(week_month_fe)
}

bg_95_logit <- temperature_deviation_bg_logit(bg_sur.s3.l, 0.95)
bg_90_logit <- temperature_deviation_bg_logit(bg_sur.s3.l, 0.90)

p_bg_logit <- ggplot() + 
  
  geom_line(data = bg_95_logit[bg_95_logit$time %in% "week", ], aes(y = estimate, x = gamma), size = 2, alpha = .35, color = "blue4") +
  geom_line(data = bg_95_logit[bg_95_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_line(data = bg_95_logit[bg_95_logit$time %in% "month", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  geom_line(data = bg_95_logit[bg_95_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 0.5, alpha = .35, color = "gray60") +
  
  geom_point(data = bg_95_logit[bg_95_logit$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, alpha = 1, color = "blue4") +
  geom_errorbar(data = bg_95_logit[bg_95_logit$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 0.4, color = "blue4", linetype = "dashed") +
  geom_point(data = bg_90_logit[bg_90_logit$time %in% "week", ], aes(y = estimate, x = gamma), size = 3, alpha = 1, color = "blue4") +
  geom_errorbar(data = bg_90_logit[bg_90_logit$time %in% "week", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = 1, size = 1, color = "blue4") +
  
  geom_point(data = bg_95_logit[bg_95_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_logit[bg_95_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_logit[bg_90_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_logit[bg_90_logit$time %in% "week_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  geom_point(data = bg_95_logit[bg_95_logit$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_logit[bg_95_logit$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_logit[bg_90_logit$time %in% "month", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_logit[bg_90_logit$time %in% "month", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 1, color = "gray60") +
  
  geom_point(data = bg_95_logit[bg_95_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_95_logit[bg_95_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = 0.4, color = "gray60", linetype = "dashed") +
  geom_point(data = bg_90_logit[bg_90_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma), size = 3, alpha = .35, color = "gray60") +
  geom_errorbar(data = bg_90_logit[bg_90_logit$time %in% "month_pars", ], aes(y = estimate, x = gamma, ymin = lower, ymax = upper), alpha = .35, size = .75, color = "gray60") +
  
  scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = bg_95_logit[bg_95_logit$time %in% "week", ]$gamma, labels=scales::percent(bg_95_logit[bg_95_logit$time %in% "week", ]$change, accuracy = 0.01))) + 
  
  ylab(expression(beta~Delta*"Temperature")) + 
  xlab(expression(gamma*" (°F)")) + 
  ylim(-0.0125, 0.065) + 
  annotate(geom = "text", x = 65:95, y = -0.0082,
           label = temp_ranges_bg,
           size = 4.15, angle = 45, fontface = 2) +
  geom_hline(aes(yintercept = 0), colour = "black", alpha = .5, linetype = "longdash") + 
  theme_bw() + 
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold", size=14),
        axis.text.x.top = element_text(angle = 90, hjust = 0),
        axis.title=element_text(size=14,face="bold"))

# Plot:

baghdad_survey_combined <- ggarrange(p_bg_logit + xlab("") + ylab("Logit"),
                                     p_bg_ols + xlab("") + ylab("LPM"),
                                  ncol = 1, nrow=2, labels  = "")

# PRODUCES FIGURE 15 IN PAPER:
annotate_figure(baghdad_survey_combined,
                top = text_grob(expression(rho), size = 18, face = "bold"),
                bottom = text_grob(expression(gamma*" (°F)"), 
                                   face = "bold", size = 18),
                left = text_grob(expression(beta~Delta*"Temperature"), 
                                 size = 18, rot = 90, face = "bold"))

###
### END.
###
